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Abstract. A coupled cell network describes interacting (coupled) individual systems (cells). 
As in networks from real applications, coupled cell networks can represent inhomogeneous networks 
where different types of cells interact with each other in different ways, which can be represented 
graphically by different symbols, or abstractly by equivalence relations. 

Various synchronous behaviors, from full synchrony to partial synchrony, can be observed for a 
given network. Patterns of synchrony, which do not depend on specific dynamics of the network, 
but only on the network structure, are associated with a special type of partition of cells, termed 
balanced equivalence relations. Algorithms in Aldis (2008) and Belykh and Hasler (2011) find the 
unique pattern of synchrony with the least clusters. In this paper, we compute the set of all possible 
patterns of synchrony and show their hierarchy structure as a complete lattice. 

We represent the network structure of a given coupled cell network by a symbolic adjacency 
matrix encoding the different coupling types. We show that balanced equivalence relations can be 
determined by a matrix computation on the adjacency matrix which forms a block structure for each 
balanced equivalence relation. This leads to a computer algorithm to search for all possible balanced 
equivalence relations. Our computer program outputs the balanced equivalence relations, quotient 
matrices, and a complete lattice for user specified coupled cell networks. Finding the balanced 
equivalence relations of any network of up to 15 nodes is tractable, but for larger networks this 
depends on the pattern of synchrony with least clusters. 
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1. Introduction. In many areas of science, interacting objects can be repre- 
sented as a network. Examples can be found in biological, chemical, physical, tech- 
nological and social systems [29, 35 . One important dynamical feature of networks 
is the possibility of synchrony, which occurs when distinct individuals exhibit identi- 
cal dynamics. Synchronization of initially distinct dynamics, such as that appearing 
in fireflies, coupled lasers and coupled chaotic systems have been extensively studied 
(see [12, 5, 39 for reviews). Many studies have investigated the role of synchrony in 
a wide range of cognitive and information processing, including recently the possi- 
ble relevance of neural synchrony in pathological brain states, such as epilepsy and 
Alzheimer's disease (reviewed in [37]). 

In this paper however we are interested in not only full synchronization, but also 
partial synchronization where a network breaks into sub-networks, called clusters, such 
that all individual systems within one cluster are perfectly synchronised. In coupled 
chaotic systems, such partial synchronization (or clustering) has been attracting great 
interest [Til H21 HOI EE] . Partial synchronies can also appear as a result of synchrony 
breaking. Suppose all individuals of a network are initially synchronized, but at some 
point lose coherence - giving synchronized sub-networks or even differently behaving 
individuals. Differentiation of (biological) cells [25] . speciation [32] . desynchronization 
of coupled oscillators [26 , or the loss of coherence in lasers can be interpreted in this 
way. We are interested in computing those potential partial synchronies which are 
solely determined by the network architecture (topology), rather that any specific 
details of the network dynamics (e.g. parameter values or function forms). 
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Mathematically such interaction networks are described as a directed graph [36J 
[41]. Nodes correspond to the individuals, and arrows (edges) denote their interactions. 
Coupled cell networks are a general formalism using a directed graph to describe 
such interacting individuals (see [34j [19j HE]). I n this settings cells correspond to 
the individuals (graph nodes), and there can be multiple types of cell, and similarly 
multiple types of arrow (graph edges). The topology of the network is described by 
an adjacency matrix, using symbolic entries for different arrow types. 

In this paper we describe how to compute all possible partial synchronies which 
are a consequence of the topology (i.e. the adjacency matrix for the graph) of a 
given coupled cell network. Such partial synchronies are represented as a partition 
of the network cells, termed a balanced equivalence relation (also referred to as a 
balanced coloring). These impose a block structure on the adjacency matrix, leading 
to a computer algorithm which determines all possible balanced equivalence relations 
using matrix computations. The set of all possible balanced equivalence relations is 
partially ordered, and forms a complete lattice (see [33). In this paper, we compute all 
possible balanced equivalence relations and their complete lattice. Existing algorithms 
in [HE] can fi n d the top lattice node, i.e. the maximal balanced equivalence relations 
from a given network topology. We use the top lattice node in order to reduce the 
search space for finding all possible balanced equivalence relations, and so speed up 
constructing the complete lattice. 

The supplementary material includes an implementation of the algorithms de- 
scribed, and a hybrid of the top lattice node algorithms from [4, 9 , using the freely 
available programming language Python ( |http:/ /www . python. org) and Numerical 
Python library (http://numpy.scipy.org) for matrix support. We hope that inter- 
ested researchers will be able to take this script and run it on networks of interest 
by entering the adjacency matrices. The code prints out balanced equivalence rela- 
tions, their quotient network adjacency matrices, and the associated lattice structure 
(as text). Additionally, provided GraphViz [15] and associated Python libraries for 
calling it are installed, figures of the network and lattice are also produced. 

Our implementation effectively finds all the balanced equivalence relations for any 
coupled cell network of up to 15 cells. For larger networks, the speed of computation 
depends on the total number of possible partitions to check, which depends on the 
clustering pattern of the maximal balanced equivalence relation, but not directly the 
size of the network. Inhomogeneous networks are most tractable, and an example of 
a 30 cell network is discussed. 

This paper is organized as follows. In Section[2] we recall some basic features of the 
coupled cell formalism, then review the basics of lattice theory. In Section [3] we show 
that finding a balanced equivalence relation is equivalent to finding a particular type of 
invariant subspace of a linear map, which is represented by the adjacency matrix of a 
given coupled cell network. We then show that an adjacency matrix, which leaves such 
subspaces invariant, has a block structure. This matrix property leads to the computer 
algorithm discussed in Section [4] which determines all possible balanced equivalence 
relations of a given network, allowing display of the corresponding lattice of the set 
of all balanced equivalence relations. Finally in Section [5] we demonstrate how the 
algorithm can be applied to several examples from the literature, with conclusions in 
Section [6] 

2. Preliminaries. 

2.1. Coupled cell network and associated coupled cell system. A coupled 
cell network describes interacting individual systems schematically by a finite directed 
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graph Q = (C,£,~c,~e)- Here C = {ci, C2 . . . , c n } is the set of nodes (cells), £ = 
{ei, e2, . . . , e m } is the arrows between them (couplings), and equivalence relations ~c 
and ~e describe different types of cells and couplings. 

Different cell types (equivalence relation ~ c ) can be labelled with symbols such 
as circles, squares, or triangles. Similarly different arrow types (~e) can be shown 
using different kinds of lines (solid, dashed, dotted). 

Each cell c is a dynamical system with variables x c in cell phase space P c , for 
simplicity a finite-dimensional real vector space R r , where r may depend on c. Cells 
of the same type must have the same phase space, for c, d G C, c ~ c d => P c = 

Each arrow e G £ connects a tail node T(e) to a head node H(e), expressed using 
maps H : £ —> C and T : £ —> C. Arrows of the same type must have matching tail 
and head cell types, ei,e2 G £, ei ~e ^2 => T!{e>i) ~c T~i{^2) an d T( e i) ~c T(e<i), 
ensuring ~e maintains similar input/output characteristics. 

For each cell c G C the input set of c is defined as 1(c) = {e £ £ : H(e) = c}, 
where e G /(c) is called an mp^£ arrow of c. 

Definition 2.1. /Tie relation ~j of input equivalence on C is defined by c ~j d 
if and only if there exists an arrow-type preserving bisection f3 : 7(c) — )► 1(d) such that 
i ~e /3(i) Vz G /(c) ; and two cells c and d are said to be input isomorphic. 

The input equivalence on C identifies equivalent cells whose input couplings are 
also equivalent. As a consequence, if c ~j d then they should have similar dynamics 
defined in a coupled cell system, which is a system of ordinary differential equations 
(ODEs) associated with a given coupled cell network. Define the total phase space 
to a given n-cell coupled cell network Q to be P = P Cl x • • • x P Cn and employ the 
coordinate system (x\, . . . ,x n ) G P, where X{ G P Ci . The system associated with the 
cell Ci has the form 



where the first variable of fi represents the internal dynamics of the cell q and the 
remaining q variables {xj 1 , . . . ,Xj } = T(I{c%)) represent coupling. In this paper, we 
employ definitions which permit multiple arrows (some subsets of indices 2k are equal) 
and self-coupling (some jk equal i). 

The function fi corresponds to the i-th component of an admissible vector field 
P = (/i 5 • • • 5 /n)? which is compatible with the network structure, and depends on a 
fixed choice of the total phase space P. It follows that different components of F are 
identical if the corresponding cells are input isomorphic, i.e., f c = fa for c ~j d. As a 
consequence, the number of distinct functions in a coupled cell system corresponds to 
the number of input equivalence classes. We now define types of coupled cell networks. 

Definition 2.2. A homogeneous network is a coupled cell network such that all 
cells are input isomorphic or identical. If a coupled cell network is not homogeneous, 
we call it an inhomogeneous network. A homogeneous network that has one equiva- 
lence class of arrows is said to be regular. The valency of a homogeneous network is 
the number of arrows into each cell. 

Table |2.1| shows various types of coupled cell networks and the corresponding 
coupled cell systems. Note that, by definition, c ~j d c ~c d holds, but the 
converse fails in general. If 1(c) = and 1(d) = 0, we say that c ~j d iff c ~c d ( see 
for example the network Q\ in Table |2?T| ) . 

2.2. Balanced equivalence relations and quotient networks. We say a 

given coupled cell system has synchrony if (at least) two cells c and d have identical 
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Table 2.1 

Coupled cell networks, and the corresponding coupled cell systems and symbolic adjacency ma- 
trices (defined in Section^. The overline indicates that some couplings from other cells to that cell 
are identical, i.e., f(xi,Xj,xj~,xi) means f(xi,Xj,xj~,xi) = f(xi,Xj,xi,xj~). 
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outputs, that is x c (t) = Xd(t) Vt G M. Synchronous network dynamics solely de- 
termined by the network structure is associated with a special type of partitions of 
cells termed balanced equivalence relations, which define a smaller network called the 
quotient network describing synchronous dynamics of the original network. 

Let ix be an equivalence relation on C, partitioning the cells into equivalence 
classes. For a given equivalence relation ex, the corresponding subspace of the total 
phase space P is defined by A M = {x <E P : c tx\ d => x c = x^}, which is called a 
polydiagonal subspace of P. 

Denote by Tq the class of admissible vector fields of a given coupled cell network Q 
with the total phase space P. A polydiagonal subspace is called a synchrony subspace 
(or balanced polydiagonal) if it is flow-invariant for every admissible vector field with 
the given network architecture. That is, 

F(A M ) C A M VF G 

Equivalently, if x(t) is a trajectory of any F G Tq, with initial condition x(0) G 
A M , then x(t) G A M for all t G R. Patterns of such robust synchrony are classified 
by a special type of equivalence relation defined in the following. 

Definition 2.3. An equivalence relation on C is balanced if for every c,d G C 
with c ex d, there exists an input isomorphism j3 such that T(i) ix T(/3(i)) for all 
i G 1(c), where the map T(e) returns the tail node of an arrow e G £. 

In particular, the existence of f3 implies c ~j d. Hence, balanced equivalence 
relations can only occur between input isomorphic cells. A necessary and sufficient 
condition for a polydiagonal subspace to be a synchrony subspace is given by: 

Theorem 2.4. An equivalence relation ix on C satisfies F(A M ) C A M for any 
admissible vector field F if and only if ix is balanced. 

Proof See [19, Theorem 4.3]. □ 

A balanced equivalence relation M on a network Q induces a unique canonical 
coupled cell network Q/^ on A M , called the quotient network. The set of cells of 
the quotient network Q/^ is defined as C M = {c : c G C}, where c denote the ix- 
equivalence class of c G C, and the set of arrows is defined as £ M = U c e<s^( c )> wnere 
S is a set of cells consisting of precisely one cell c from each ix-equivalence class. 

For example, Figure |2.1| shows all balanced equivalence relations of the inho- 
mogeneous network C?3 in Table |2.1| and the corresponding quotient networks. Any 
dynamics on the quotient lifts to a synchronous dynamic on the original network Q. 
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txii= (13)(245) m 2 = (13)(24)(5) m 3 = (1)(25)(3)(4) ^ 4 = (1)(2)(3)(4)(5) 



Fig. 2.1. All balanced equivalence relations of the inhomogeneous network Q3 (see Table \2.1y and 
the corresponding quotient networks. Note that the balanced equivalence relation 1x14= (1)(2)(3)(4)(5) 
is trivial and this gives the original 5- cell network. 



2.3. Lattice theory. All possible partial synchronies (balanced equivalence re- 
lations) have a hierarchy structure represented as a complete lattice. We recall some 
basic facts about lattice theory using balanced equivalence relations as an example 
for some concepts. See [14] for concepts in general and for more details. 
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The set of all balanced equivalence relations has a partially ordered structure, 
using the relation of refinement. Let ix^ and Mj be balanced equivalence relations on 
the set C. Recall that Mi refines Mj, denoted by ix^lXj, if and only if c txi^ d => c Mj d 
where c, d G C. That is, [c]i C [c]j where [c]k is the IX/- -equivalence class. The set of 
all balanced equivalence relations of a (locally finite) network form a complete lattice 
in general [33j H] . 

A complete lattice has a top (maximal) element, denoted T, and bottom (minimal) 
element, denoted _L For example, the top element of the complete lattice of balanced 
equivalence relations for any n-cell homogeneous network is trivial and given by Mj= 
(12- — n) (i.e., all cells are synchronous). Aldis [4] and Belykh and Hasler [9] give 
algorithms to find a nontrivial maximal balanced equivalence relation (top). For any 
n-cell coupled cell network, the bottom element is !X^ = (1)(2) • • • (n) (i.e., all cells 
are distinct). 

The structure of a lattice can be visualised by a diagram. Let DX^, Mj and ix& be 
distinct balanced equivalence relations. We say tx^ is covered by iXj, denoted M^<iXj, 
if and only if DXLj^lXj and ix^ix^iXj holds for no IX&. In a diagram, circles represent 
elements of the ordered set, and two elements ix^, iXj are connected by a straight line 
if and only if one covers the other: if is covered by IXj, then the circle representing 
ix j is higher than the circle representing ix$. The rank of an equivalence relation is 



the number of its equivalence classes (see [23]). Figure 2.2 shows the complete lattice 



of the partially ordered set of all 4 balanced equivalence relations of the coupled cell 



network Qs, which are listed in Figure 2.1 



Rank 1 




Fig. 2.2. The complete lattice of the partially ordered set of all possible balan ced equivalence 
relations of the coupled cell network Q3 along with their ranks (see also Figure \2A^ . For example, 
the balanced equivalence relation 1X3= (1)(25)(3)(4) has 4 equivalence classes, hence its rank is 4. 
Since this network is inhomogeneous , the top element is not of rank 1, which requires all cells are 
identical. 



3. Matrix computation for balanced equivalence relations. We aim to 

determine robust patterns of synchrony of a coupled cell network solely from the 
network structure, which is described by the corresponding symbolic adjacency matrix 
defined as follows. 

Definition 3.1. Let Q = (C,£ ^c^e) be an n-cell coupled cell network with I 
cell-types and m arrow-types with [ci]c, • • • , [q]cs the ~c -equivalence classes for cells 
and [ei]#, . . . , [e m ]# ; the ~e -equivalence classes for arrows. We define the symbolic 
adjacency matrix ofQ to be the nxn matrix A = (a^). The (i, j) -entry corresponds to 
the number of arrows of types [ei]#, . . . , [e m ]# from cell j to cell i, represented by the 
sum Y^k=i Pk^k, where e^ is the type of arrow corresponding to the [e^Js- equivalence 
class and (3k is the number of arrows of type e^. 
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Example coupled cell networks and their corresponding symbolic adjacency ma- 
trices are shown in Table I2TT1 

Our main result determines all possible balanced equivalence relations combina- 
torially using matrix manipulations. In Lemma [3~8t we show that the necessary and 
sufficient condition for a synchrony subspace imposes a matrix property defined as 
follows. 

Definition 3.2. Let B = (6^-) be a p x q symbolic matrix. We say B is a 
homogenous block matrix if the sum Sj=i ^ij ^ s identical for all rows i = 1, . . . ,p. 

The polydiagonal A M is defined by a given equivalence relation \x\ which deter- 
mines a unique partition of cells. We use normal form cycle notation which is obtained 
by writing the cell numbers 1, . . . ,n in increasing order in each cycle, starting with 
the 1-cycle, then the 2-cycles, and so on in increasing order of length. For example, 
the following polydiagonal subspace A M = {(xi, #3, #4, #5, #6)1^2 = £4,^3 = ^5 = 
xq } corresponds to the equivalence relation 1x1= (1) (24) (356) and can be written as 
1x1= [1 1 2 1 3 1 ]. Let A = (aij)j i,j = 1, . . . , 6 be the adjacency matrix of a 6-cell coupled 
cell network. For the above normal form cycle notation of the equivalence relation ixj, 
we arrange the columns and rows of the adjacency matrix of the network accordingly 
as follows: 
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This is a block matrix, with 3x3 blocks, where 3 is the number of equivalence classes. 

More generally, let txi= [l ai 2 a2 • • • n an ] and k = J^ILi a i be the number of equiv- 
alence classes of M which determines the polydiagonal A M . Interchanging the adja- 
cency matrix rows and columns to match the permutation normal form gives a block 
matrix with k x k blocks. 



Our main result (Theorem 3.9) states that a polydiagonal subspace A M is a 
synchrony subspace if and only if each block of the (reordered) adjacency matrix A is 
a homogeneous block matrix. 

Alternatively, our result can be obtained by defining one (integer entry) adjacency 
matrix per arrow type (as in [3] for regular networks with one arrow type), and finding 
the intersection of balanced equivalence relations for the arrow-specific adjacency 
matrices. 



3.1. Linear admissible vector fields. We represent the matrix form of Q- 
admissible linear vector fields on an n-cell coupled cell network Q. 

Definition 3.3. Let S = (5^) and T = (tij) be symbolic matrices with the 



same size. Suppose Si 



j/ for some (i,j)-th and (i',j')-th entries of S. If the 



corresponding entries in T satisfy Uj = tvy for all such indices (i,j) and then 
we denote this relation by S ~ T. 

Proposition 3.4. Let A = (a^) be the n x n symbolic adjacency matrix of an 
n-cell coupled cell network Q. The n x n symbolic matrix J = (Jij) of Q -admissible 
linear vector fields on Q has the form: 



J = D + A, 



(3.1) 
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where D = (dij) is an n x n diagonal matrix with da = dkk when i ~j k for i,k G C, 
and A = (dij) is an n x n matrix such that A ~ A. 

Proof. Let t\ be the dimension of internal dynamics of the z-th cell. Let J be the 
n' x n' matrix form of the ^-admissible linear vector field of the network Q with total 

n 

phase space P' = R n where n' = ri . 

i=l 

J can be described as a block matrix J = (Jij), i, j = 1, . . . , n, and each block of 
Jij is a hi x kj matrix with real entries. We can decompose J as 

J = D + (-D + J). 

Each block of the diagonal matrix Da satisfies Da = Dkk when i ~j k for i,k G C. 
By representing each block Da with a symbol da, we obtain annxn symbolic diagonal 
matrix D. Similarly, by representing each block of —D + J with a symbol a^, 
we obtain an n x n matrix A, which satisfies A ~ A by the admissibility. Therefore, 
the n x n symbolic matrix J = (J^) of ^-admissible linear vector fields on an n-cell 
coupled cell network Q has the form J = D + A. □ 

An equivalence relation ixi is balanced if and only if the ^-admissible linear vector 
field J satisfies J(A M ) C A M [19, Theorem 4.3]. This fact leads to a generalization 
of Proposition 4.1. in [23], which consists of only one type of cell and one type of 
coupling. 

Proposition 3.5. Let Q he a coupled cell network associated with the admissible 
vector field F. Let A be the adjacency matrix of Q. M is balanced if and only if 
A(A M ) C A m where A M is a synchrony subspace associated with dxl 

Proof Using Equation |3.1[ the result follows immediately: 

J (A[xi) C A M 
^{D + i}(A M ) C A M 
o{D(A N )+I(A N )}CA N 

x since ixi refines ~/ , and therefore D(A M ) C A M 
<^> A(A M ) C A M since A ~ A. 

□ 

Now, we aim to determine all possible invariant polydiagonals under a given adja- 
cency matrix using a projection map. Since the determination of synchrony subspaces 
does not depend on the size of internal dynamics of the cells, without loss of gener- 
ality, we assume the total phase space P = W 1 for an n-cell coupled cell network for 
the remaining arguments. 

3.2. Projection onto a polydiagonal. We construct a projection map on a 
given polydiagonal as follows. 

Let A M C R n be a polydiagonal subspace of R n , and A^ denote its complement. 
We define the projection map P M of R n on A N along A^ by Im(P M ) = A M and 
Kvri I ',)-/■,'... 

Let A M be the polydiagonal determined by a given equivalence relation ix. We use 
normal form cycle notation for an equivalence relation. For example, the equivalence 
relation M corresponding to a polydiagonal subspace A M = {(#i,#2>#3)|#2 = ^3} is 
written as a product of disjoint cycles 1x1= (1)(23) in normal form, and also written 
as m= [l^ 1 ]. 
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Now we define a map n which maps each element to the first element of the cycle 
that they belong to when written in normal form. For example, the elements in the 
above product of disjoint cycles are mapped to 1 — >• 1, 2 — >• 2, and 3 — » 2 by a map tt. 
We define the corresponding projection matrix P M = (p^-) on A M , which is written 
in normal form of a partition using the map tt as follows: 

with all other entries being 0. 

With the elements of P M defined in this way, the projection matrix has a block 
diagonal form: 

f P 1 • • • \ 
P 2 ••• 

V ••• P k J 

where k is the number of disjoint cycles and P^, i = 1, • • • , k is a t\ x t{ square 
projection matrix on A = {(^i, . . . ,x ti )\xi — • • • = x ti } with rank(P^) = 1 Vz and 
off-diagonal blocks are zero matrices. 

Lemma 3.6. Let P M and A be linear mappings ofW 1 and let W 1 = A M ® A^. 
A M is A-invariant if and only if P M AP M = AP M; where P M is the projection on A M 
along A^ and A is the adjacency matrix of a given coupled cell network. 
Proof This result is well known, see for example [27] and [23]. □ 
For the rest of arguments, we arrange the columns and rows of the adjacency 
matrix of the network according to the normal form of a given equivalence relation 

1X1. 

Proposition 3.7. A M is a synchrony subspace of a coupled cell network Q if 
and only if P^AP^ = AP M; where P M is the projection on A M along A^ and A is 
the adjacency matrix of Q . 

Proof If A M is a synchrony subspace, then the corresponding equivalence relation 
ixi is balanced. Hence, by Proposition |3.5[ A M is A-invariant and from Lemma [3.6[ 
the corresponding projection map P M to A M satisfies P M AP M = AP M . Conversely, if 
P^AP^ = AP M , then the corresponding subspace A M is A-invariant by Lemma |3l)] 
and Proposition |3.5[ thus the corresponding polydiagonal is balanced. □ 



3.3. Block structure of an adjacency matrix. We show that the necessary 



and sufficient condition for a synchrony subspace in Proposition 3.7 imposes a block 
structure on the adjacency matrix. 

Lemma 3.8. Let A be the nx n adjacency matrix of a given coupled cell network 
Q . Suppose a synchrony subspace A M is defined by a partition [l ai 2 a2 • • • n an ] and the 
corresponding projection matrix P M is a block diagonal matrix whose diagonal blocks 
Pi, i = 1, . . . , k, where k = J^ILi a i are projection matrices on diagonal subspaces A 
in the corresponding dimensions. Then P^AP^ = AP^ if and only if corresponding 
blocks of A to P M satisfy the following condition: 

• Let A st , where s, t = 1, . . . , k be blocks of A. For all blocks A st , the sum of 
each row is identical. 
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/ Pi 



Proof. Since P M = 



AP^ — 



P^AP^ — 



V o 



( Pi 



V o 



we obtain: 



Pk J 
A lk \ ( Pi 



\ Ak\ ■ ■ ■ Akk J 
( AuPj A 12 P 2 •• 
A 2 \P\ A22P2 ■■ 



V ••• 

A lk P k ^ 

A 2k P k 



\ 
Pk J 



\ A kl P x A k2 P 2 







• • • A kk P k J 

( AnPi A 12 P 2 
A 2 \P\ A 22 P 2 



\ A kl P x A k2 P 2 

( PiAnPi PiA 12 P 2 • • • P x A lk P k \ 
P2A 21 P 1 P 2 A 22 P 2 ••• P 2 A 2k P k 



\ P k A kl P x P k A k2 P 2 



PkA kk P k J 



A lk P k \ 
A 2k P k 



A kk P k ) 



Hence, P^AP^ = AP M P s A st P t = A st P t , for all s, t = 1, . . . , fc. 

Let A st = (a st )ij be an arbitrary I x m block matrix. Then P t is a m x m square 
matrix and P s is a Z x I square matrix. Since P t and P s are projection matrices onto 
m-dimensional and /-dimensional diagonals, respectively, 



A st Pt 



( 1 
1 







P,A, t P, 



< 



\ off + ' 



/ 1 



V 1 

/ <i + 

«?i + 

V < + 



lm 

\ 




+ a 



lm 



V 1 



/ 



<m 



\ 




0/ 



+ a 



lm 
1 st 
l 2m 



J V ag + --- + a 

••• \ 

••• 

••• J 



St 

lm 







/ 



Thus P s A st P t = A s tPt if and only if Y^jLi a tj * s identical for alH = 1, . . . , I. 
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Therefore, P s A st P t = A st P t for all s, t = 1, . . . , k if and only if, for all blocks A st , 
the sum of each row is identical. □ 

We call a block of this form a homogeneous block matrix. It now follows immedi- 
ately that: 

Theorem 3.9. A poly diagonal subspace A M is a synchrony subspace if and only 
if each block of the adjacency matrix A, which corresponds to a block of P M; is a 
homogeneous block matrix. 

Proof. Each block of A is a homogeneous block matrix 
P M AP M = AP^ (by Lemma [3?8| 
^==> A M is a synchrony subspace (by Proposition |3.7[ ) . □ 

Example 3.1. The projection mapping on A M = {(xi, x 2l ^3)^2 = ^3} has the 
form: 

















1 









1 






Let A be the adjacency matrix of a 3-cell coupled cell network of the form: 





( an 


«12 


«13 




&21 


«22 


«23 




\ «31 


«32 


«33 



Then, 



P^AP^ 




an 


d!2 


«13 ) 




«21 


«22 


«23 


1 


«31 


«32 


«33 




an 


«12 " 


f «13 





&21 


&22 - 


f «23 





«31 


&32 - 


f «33 








an &12 

&21 &22 
&21 &22 

Therefore, 

P M AP M AP^ <=> a 2 i a 3 i and a 2 2 + «23 «32 + «33- 
Hence, A M zs a synchrony subspace if and only if A has the following block structure: 





/ an 


ai2 


«13 


H 




«22 


«23 




\ «21 


«32 


«33 



w/iere a 2 2 + «23 = «32 + «33- ^ 
Next we derive the adjacency matrix of the quotient network which is defined as 

the adjacency matrix A restricted on a synchrony subspace A M . 

Corollary 3.10. Let A M be a synchrony subspace defined by a partition [\ OL1 2° 1 ' 2 • • 
with k = ct\ + • • • + a n equivalence classes, and P M be the corresponding block projec- 
tion matrix. Let A st for s,t = 1, . . . , k corresponding to the blocks of P M be blocks of 
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an n x n adjacency matrix A. If blocks A st = (a st )ij are homogeneous block matrices 
such that: 



( M 



A± ai 



A 



iQi+i) 



^l(ai + -+Q»-i 



±11 



A lk \ 



i21 



i2ai 



i 2(ai+Qj 2 ) 



i 2(ai + -+a Tt -i + l) 



A 2k 



Akk J 



^fc(ai + l) 



A 



fc(ai+a 2 ) 



^(aii-.+an-i+l) 



tften £fte quotient network corresponding to ix /ias a kxk adjacency matrix A|a m , 
denoted by A^, of the form: 



En lOiH harj 
.7 = 1 a lj 



1 + 1) 



2(ai + ---+a n _i + l) 



Sz=i a ij 



V 1 o 2ai 



En 2(< 



En 
2=1 



,2fc 



n fe(aiH h^r, 

j=l a lj 



El 
i=i a ij 



Proof. Let {i> i, . . . , v ai , . . . , t^} be a basis of a synchrony subspace. Each 

basis element corresponds to a conjugacy class of the partition ix and is an n x 1 
vector. Therefore, the basis elements have the following forms: 



vi = [1,0,... ,0]* 
v ai = [0 1 _^ I 0,l,0,...,0] t 

Oi\—l 

v ai+1 = [0 1 __0,^,0,...,0] t 

on 2 



v fc = [0,...,0,l,...,l] 



Since each block A s£ is a homogeneous block matrix, i.e., the sum of each row is 
identical, we can express the image of each basis element using a linear combination 
of a basis with the sum of the first row of each A st being used as a coefficient such as: 



COMPUTATION OF BALANCED EQUIVALENCE RELATIONS 



13 



Avi 



Avo 



i i 

V 1 a 22 



1 1 



E 



fr9 



En 
.7 = 1 



X lj 



n 



E 



CLljVk 



Therefore a k x k matrix A M , which is the adjacency matrix A restricted on A M , 
is written as 



2^7 = 1 a i7 



V 1 a 2ai 

^.7 = 1 a l.7 



e.Li « 



En l(o:iH \-a 
.7 = 1 a l.7 



i + l) 



2(ai + -"+a„_ 1 + l) 
^ 



n /c(aiH ha^-i + l) 

J=l a lj 



En 
If 



2fc 
1 a lj 



and this is the adjacency matrix of the quotient network corresponding to IX. □ 
4. Computer algorithms. 

4.1. Balanced equivalence relations. The above combinatorial properties of 
adjacency matrices leads to a computer algorithm which determines all balanced 
equivalence relations and adjacency matrices A M of associated quotient networks 
for a given coupled network Q. 

We enumerate all possible equivalence relations x of n-cells, and test which are 
balanced. If the top lattice node is known in advance, e.g. using the algorithm in 
[4] or [9], only equivalence relations which refine this top node need be tested. To 
test if each ix is balanced, we construct an n x k matrix, where k is the number 
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of equivalence classes of ixi, from the n x n adjacency matrix of Q. All rows in each 
equivalence class are identical for a balanced equivalence relation. Finally, for balanced 
equivalence relations, we can construct adjacency matrices A M of the corresponding 
quotient networks Q/^. 

Step 1: For a given n-cell coupled cell network Q = (C,£,~c,~e), we express the 
corresponding n x n adjacency matrix A as 

A=[d---C n ] 

where C{ G R nxl , i = 1, . . . , n are column vectors. Let C p denote the M-equivalence 
classes on C where p = 1, . . . , k. For example, if C = {{1, 3, 5}, {2}, {4}} then C\ = 
{1,3,5}, C 2 = {2}, and C 3 = {4}. Note that J2p=i \C P \ = n. Let C pl be the first 
element of each equivalence class. We assume that C\\ < C 2 i < • • • < C/ci, where 
these cell numbers are used as indices jbr row vectors in Step 3. 
We generate a new n x k matrix A M with columns 



jec p 



for all possible equivalence relations m. 



Let R 



plxk 



where 



for p ■ 



denote the row vectors of this new n x k 



matrix Therefore, 



/ Ri 



'[X]! 



Step 2: Now we determine which equivalence relations ix are balanced. An equiva- 
lence relation ix on C is balanced if and only if for all p = 1, . . . , k we have: 



Rt 



■[Xll 



Rt> 



V/, m G C p 



(4.1) 



Step 3: If the above condition 4.1 is satisfied, the k x k adjacency matrix of the 



quotient network A M corresponding to a balanced equivalence relation ix is given by: 

/ RtXi! 



where G 



l lxk , i = 1, . . . , k are representative row vectors in each equivalence 



class. 



Example 4.1. Consider the homogeneous network Q§ in Table 2.1 with the cor- 



responding adjacency matrix shown in Figure J^l 



We determine if the equivalence relation (i.e., partition of cells) ix= (135) (24) 
is balanced or not by the matrix computation. There are two equivalence classes 
C\ = {1,3,5} and C2 = {2,4}. We generate a new 5x2 matrix by adding 
columns 1, 3 ; 5 and 2, 4 such that 



( 

2ei + e 2 
ei + e 2 
\ 2ei + e 2 



2ei + e 2 \ 



ei 
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/ ei+e 2 ei \ 

ei ei + e 2 

e 2 ei ei 

e 1 +e 2 ei 

\ e 2 ei e 1 J 



Fig. 4.1. Homogeneous network Q§ in Table pO] with the corresponding symbolic adjacency 
matrix. 



The equivalence relation x is balanced if and only if 
[0 2ei + e 2 ] = [2ei + e 2 0] = [2e x + e 2 0] ana 7 [ei ei + e 2 ] = [ei + e 2 ei]. 
However, this does not hold. Thus the equivalence relation ix= (135) (24) zs no£ 6a/- 

On £/ie o£/ier /mna 7 , Ze^ ixi= (124) (3) (5). There are three equivalence classes C\ — 



{1,2,4}, C 2 = {3} and C 3 
columns 1, 2, 4 such that 



{5}. We generate a new 5x3 matrix A M 6 2/ adding 



( 2ei + e 2 
2ei + e 2 
e 2 

2ei + e 2 
V ^ 






ei 




\ 




ei / 



77ie equivalence relation ix zs balanced if and only if 
[2ei + e 2 0] = [2e x + e 2 0] = [2e x 



e 2 0]. 



This is satisfied. Thus the equivalence relation ix= (124) (3) (5) zs balanced. As a 
result, the quotient network corresponding to the balanced equivalence relation 
1x1= (124)(3)(5) and the associated 3x3 adjacency matrix A M are given z'n Figure 4-2, 



QU 



2ei + e 2 
e 2 ei ei 
e 2 ei ei 



Fig. 4.2. T/ie quotient network G/\x\ corresponding to a balanced 
(124) (3) (5) and the associated 3x3 adjacency matrix A^. 




ence relation 1X1= 





The algorithm as shown above describes the graph Q with a single adjacency ma- 
trix A containing symbolic entries for the different arrow types. We now discuss an 
alternative representation using separate integer matrices for each arrow type, which 
for most programming languages is more practical to implement. The following defi- 
nition is a variation of Definition 5.2 in [2] for homogeneous networks. 

Definition 4.1. Let Q = (C,£,~c,~e) be an n-cell coupled cell network with I 
cell-types and m arrow-types with [ci]c, • • , [q]c> the ~ c - equivalence classes for cells 
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and [ei ]#,... , [e m ]# ; the ~e -equivalence classes for arrows. We define the adjacency 
matrix of Q with respect to [ey^E, for k = 1, . . . , m to be the nxn matrix M^g^) • The 
(i, j) -entry corresponds to the number of arrows of types [ej^E from cell j to cell i. 
Notice by construction we have: 



k=i 



Therefore the above algorithm procedure can now be applied to each of the m 
arrow type specific matrices individually, and if it holds for all of them, it also holds 
for A as well. We now repeat Example |4.1| to demonstrate this. 

Example 4.2. For the coupled cell network Q , two adjacency matrices M^g^ = 
(mjj) (solid) and M^g^) = (j^fj) (dashed) for two different arrow types are defined as 
in Figure 




M, 



( 


1 





1 


o \ 




( 


1 








o \ 


1 








1 
















1 











1 





1 




1 














1 


1 











1 














V o 





1 





1 J 




V 1 











o / 



Fig. 4.3. Homogeneous network Q§ in Table 2.1 with two adjacency matrices M^g^ = (m|.) 



(solid) and M(g 2 ) = (dashed) for two different arrow types. 



Using these arrow type specific adjacency matrices, we determine if the equivalence 
relation ix= (135)(24) is balanced. There are two equivalence classes C\ = {1,3,5} 
and C2 = {2,4}. We generate new 5x2 matrices by adding vectors in columns 1, 3, 
5 and 2, 4: 



M< 



( 2 \ 

1 1 

2 
1 1 

\2 0/ 



( 1 \ 

1 

1 
1 

V 1 / 



The equivalence relation x is balanced if and only if for each arrow specific 5x2 
matrix rows 1, 3 ; and 5 are equal, and rows 2 are 4 equal. However, this does not 
hold. Thus the equivalence relation ix= (135) (24) is not balanced. 

On the other hand, let 1x1= (124) (3) (5). There are three equivalence classes C\ — 
{1,2,4}, C2 = {3} and C3 = {5}. We generate new 5x3 matrices by adding vectors 
in columns 1, 2, 4: 



/ 2 

2 

2 

V 



\ 


1 



1 / 



Mi 



( 1 
1 
1 
1 

V 1 



\ 




/ 
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The equivalence relation dxi is balanced if and only if for each arrow specific 5x2 
matrix rows 1, 2, and 4 are equal. This is satisfied. Thus the equivalence relation 
ixi= (124) (3) (5) is balanced. As a result, the quotient network £//tx corresponding to 
a balanced equivalence relation dxi= (124) (3) (5) and the associated 3x3 arrow type 
specific adjacency matrices are given in Figure 

At><\ Q I x 




(5,2) 






Fig. 4.4. The quotient network Q/^ corresponding to a balanced equivalence relation \xi= 
(124) (3) (5) and the associated 3x3 arrow type specific adjacency matrices. 



4.2. Lattice of balanced equivalence relations. Using the above computer 
algorithm, we can determine all balanced equivalence relations and corresponding 
quotient networks for a given coupled cell network. Now, using the refinement relation, 
we construct a complete lattice of balanced equivalence relations for a given coupled 
cell network. 

Let p be the total number of balanced equivalence relations of a given coupled 
cell network. We aim to compute a p x p adjacency matrix L = (Z^) for the lattice 
with entries 1 where ixi^ is covered by txjj, and otherwise. 

Step 1: Without loss of generality, order the p balanced equivalence relations by 
increasing rank (number of equivalence classes). This ensures that the top element is 
first and the bottom element last, and that the matrix L will be lower triangular. 
Step 2: Construct p x p matrix B — (pij) with entries 1 where DX^Mj (ix^ refines 
\x\j) and otherwise. This is almost the desired adjacency matrix, but it includes 
extra edges since refinement is not as strict as covering. 

Step 3: Calculate p x p matrix T = (tij) = B 2 . Non-zero entries Uj indicate nodes i 
and j are connected by a path of length two via some intermediate third node fc, thus 
txi^M/c^ixij, meaning ixi^ is not covered by Mj. We can assume k is distinct from i 
and j since the diagonal entries of B are zero. 

Step 4: Construct p x p matrix L = (Z^) using 1^ = 1 if b^ = 1 and t^ — 0, and 
kj = otherwise. 

For larger lattices computing the full matrix T in step 3 is increasingly time 
consuming, and only a fraction of the values are needed in step 4. For computational 
efficiency, only where bij = 1 do we need to check if Uj = 0. We do this by considering 
the existence of a two step path between lattice nodes i and j via node k (i.e. b^ = 1 
and bkj = 1), and as a further optimization only those nodes k with rank(z) < 
rank(/c) < rank(j) need be considered. 

The matrix L is the adjacency matrix for the p-node lattice, and defines the set 
of edges. Lattices are by convention drawn as diagrams with an up/down orientation 
with the top lattice element higher than the bottom lattice element. Additionally we 
require lattices nodes of the same rank (number of equivalence classes) to be shown 
at the same height. 



Example 4.3. Consider the five-cell homogeneous network in Table 2.1. This 
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network has 5 balanced equivalence relations as shown in Figure \475[ in rank order. 
Qs Balanced equivalence relations 



(12345) 

(124)(35) 

(124)(3)(5) 

(1)(2)(35)(4) 

(1)(2)(3)(4)(5) 




Fig. 4.5. Five-cell homogeneous network Q§ in Table \2A\ and all possible 5 balanced equivalence 
relations. 

We construct the 5x5 matrix B, which represents refinement relations between 
the 5 balanced equivalence relations. Then T = B 2 is a simple matrix multiplication, 
which is used to remove unwanted edges from matrix B to give L. Table \4-l] shows 
these three matrices. Matrix L is used as the adjacency matrix when drawing the 
lattice, with vertical positions dictated by the lattice node ranks (Figure 4-6 right). 



B 



T = B 2 



( \ 

1 

110 

110 

\ 1 1 1 1 J 



( \ 



1 
1 

\ 3 2 / 



/ \ 

1 

10 

10 

\ 1 1 ) 



ED 



Table 4.1 

Three 5x5 matrices, B , T , and L defined in the lattice algorithm for the network Q§ in Table 



Covering relations 


Lattice of b 




Rank 1 


(124) (35) < (12345) 


i 


(124)(3)(5) < (124)(35) 


Rank 2 


(1)(2)(35)(4) < (124)(35) 


Rank 3 


(1)(2)(3)(4)(5) < (124)(3)(5) 


i 

Rank 4 


(1)(2)(3)(4)(5) < (1)(2)(35)(4) 


! 




Rank 5 




Fig. 4.6. The covering relations and the lattice of balanced equivalence relations of the five-cell 
homogeneous network Q§ in Table \2A\ 



Computing all the balanced equivalence relations of a network size n scales with 
the number of equivalence relations, given by the Bell number, and is thus combina- 
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torial with the network size n. On a recent computer (2008 Apple Mac Pro) using 
a single CPU using the brute force approach with a single edge type, small networks 
take less than a second to compute, 10 nodes about 20 seconds, 11 nodes about 2 
minutes, 12 nodes about 15 minutes, 13 nodes under 2 hours, and 14 nodes about 12 
hours (with variation depending on the network topology). 

We have also implemented the algorithm in [9], and generalized this to consider 
multiple arrow types. The result is equivalent to a single-phase simplification of the 
algorithm in [4 , but less complicated to implement (see appendix). Where we have 
computed the full list of balanced coloring and identified the unique minimal balanced 
coloring, the results agree. As described above, this offers a shortcut when computing 
all the balanced equivalence relations, although for highly symmetric networks this 
optimization has limited benefit - and for regular networks offers no improvement. 
However, the time saving can be dramatic especially for random networks. As an 
example, a bidirectionally coupled chain (where the end nodes do not have self cou- 
pling) of up to 20 nodes takes under a second, 30 nodes takes about 20 seconds, and 
40 nodes about 10 minutes. 

The time (and the memory requirements) needed to compute the lattice scales 
quadratically with the number of lattice nodes. In the worst case of a fully connected 
network all equivalence relations are balanced, giving the largest possible lattice and 
the longest compute time, taking about a second for the 877 node lattice (n = 7), 20 
seconds for the 4140 node lattice (n = 8), and 10 minutes for the 21147 node lattice 
(n = 9). 

In short, in our current algorithm implementation computing the balanced color- 
ings of regular networks more than 15 nodes is impractical, although inhomogeneous 
networks are much easier to deal with. Additionally, the computations could in prin- 
ciple be run in parallel across multiple CPU cores, giving a potential linear speed 
up. 

5. Examples. The lattice of partial synchronies computed by the algorithm 
shown tells us about the existence of all possible partial synchronies determined by 
the given network structure. In this section, we select example topics from synchro- 
nized chaos and coupled neuron models, and demonstrate how a symbolic adjacency 
matrix can be defined for each example problem, and construct a complete lattice 
of all possible partial synchronies derived from the given networks structure. Some 
dynamical properties such as the stability of possible partial synchronies depends on 
the specific form of the given vector field. We demonstrate the numerical analysis of 
the stability of partial synchronies for the topic of synchronized chaos. 

5.1. Coupled identical Rossler systems. We consider a bidirectional ring of 
six diffusively coupled Rossler systems = (xi,yi,Zi) Gl 3 for i = 1, 2, . . . , 6: 

±i = -{yi + Zi) + e(^_i - 2xi + x i+1 ), 

Vi = Xi + ayi, 

%i — b + (xi — c)zi, 

with periodic boundary conditions Xq = Xq and x^ = X\. 

Since this is a regular network, the adjacency matrix consists of no n- negative 
integers and the admissible vector field is defined by a single map f , which is realized 
by the above defined system. Table [5TT] shows the adjacency matrix and the associated 
coupled cell system. The complete lattice of balanced equivalence relations of this 
network is given in Figure [5TT] 
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adjacency matrix coupled cell system 



A = 



( 


i 











1 \ 


Ul 


= f(ui,u 6 ,u 2 ) 


1 





1 











U2 


= f(u 2 ,ui,u 3 ) 





1 





1 








U3 


= f(u 3 ,u 2 ,u 4 ) 








1 





1 





u 4 


= f(u 4 ,u 3 ,u 5 ) 











1 





1 




= f(u 5 ,u 4 ,u 6 ) 


V 1 











1 


o i 


u 6 


= f(u 6 ,u 5 ,ui) 



Table 5.1 

Adjacency matrix which represents interactions among Rossler systems and the associated cou- 
pled cell system. The overline indicates that influence from coupling cells to that cell are identical, 
i.e., f(xi,Xj,x k ) means f(x i: x j: x k ) = f(xi,x k ,Xj). 




Fig. 5.1. Lattice of balanced equivalence relations of bidirectional ring of six diffusively coupled 
Rossler systems. Some partial synchronies are symmetrically related. For example the balanced 
equivalence relations (1245) (36), (1346) (25), and (14) (2356) are permutation symmetric and they 
give the same pattern of partial synchrony (a, a, f3, a, a, (3) in a bidirectional ring. 



In the numerical analysis, we take parameter values a = 0.2, b = 0.2, c — 5.7, 
and vary the coupling parameter e for the stability analysis of synchrony subspaces 
of this specific vector field. The full synchrony subspace A = {ui = • • • = uq} is 
globally stable from e ~ 0.2 to e « 1.0, and attracts all trajectories starting from 
randomly chosen initial conditions. With the loss of stability of full synchronization 
above e « 1.0, the synchrony subspace A M = {ui = u 3 = 115, u 2 = u 4 = ue} becomes 
globally stable, and attracts all trajectories. Figure 5.2 (a) shows the behaviour of 
x\ — x 2 , which is the difference between the first internal variables of the Rossler 
systems ui and u 2 , when changing the parameter e. This figure shows the gain and 
loss of stability of the full synchrony subspace. Figure [5T2] (b) shows the behaviour of 
x\ — x 3 , where the corresponding variables Ui, u 3 remain synchronized above e « 1.0. 
Only the behaviors of x\ — x 2 and x\ — X3 are illustrated in Figure pT2l qualitatively 
the same behaviors are observed for the other pairwise comparisons. 

Rossler systems with x-component coupling, as shown in this example, are known 
to exhibit a phenomena called short wavelength bifurcations [20] in which the syn- 
chronous chaotic state loses its stability with an increase of coupling strength. The 
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desynchronization behaviour of three diffusively coupled Rossler systems with Neu- 
mann boundary conditions (a bidirectional chain with self-coupling of the end cells) 
was analysed in [TT]. They showed that the singularity of the individual Rossler sys- 
tems and the use of x-coupling impose the existence of equilibria that lie outside of 
the fully synchronous subspace for any coupling strength, and proposed a direct link 
to the mechanism of desynchronization. 

A lattice-like hierarchy of synchrony subspaces of diffusively coupled identical 
systems was also discussed in [TT]. For the chain network associated with Neumann 
boundary conditions, they hypothesized a clustering type hierarchy structure based 
on the number of nodes n and its divisors, which we have verified for up to n = 15 
(see supplementary material). A related result for a linear chain with feedback in [33] 
gives an explicit lattice construction. 

5.2. Coupled Lorenz systems with heterogeneous coupling. In the next 
example, we demonstrate the symbolic adjacency matrix can be interpreted not only 
as a network structure, but also as different coupling strengths of identical individual 
systems. Consider cluster synchronization in an ensemble of five globally coupled 
Lorenz systems u^ = (x^ Zj) G R 3 for i = 1, • • • ,5 with heterogeneous coupling: 




Vi =Xi(p- Zi) -yi 

Zi — XiUi @Zi 



where a = 10, p = 28, /3 = 8/3, Ni = Y^j=i an d the coupling matrix G = (g%j) is 
defined as 



G 



( 1 1 2 4 \ 
10 12 4 
110 15 
1 2 3 2 
\ 3 2 1 2 J 



We may regard an integer value g^ to be a weight from Lorenz system j to 
Lorenz system z, and different non-zero integer values to be different weights from 
the corresponding systems. Note that Ni = 8 for all i = 1, . . . , 5. Then symbolically, 
the above coupling matrix G can be represented with five different symbols, a, 6, c, 
and e. Table [5^2] shows the symbolic adjacency matrix and the associated coupled 
cell system. 

Figure |5.3| shows five globally coupled Lorenz systems with different weighted 
arrows represented in different colors and the associated lattice of balanced equivalence 
relations of the network. If all weights are identical, the corresponding lattice of partial 
synchronies is the same as the partition lattice of 5 elements with the Bell number 
B$ = 52 lattice points. However, in this example with non-identical weights, there 
is only one non-trivial balanced equivalence relation given by (12) (3) (4) (5), which is 
found to be unstable by numerical analysis. 

We note the intriguing phenomenon termed bubbling [6 is related to the stability 
of synchrony subspaces. When the dynamics on the synchrony subspace is a chaotic 
attractor, small perturbations along the transverse direction of the synchrony subspace 
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Fig. 5.2. The stability analysis of synchrony subspaces of a ring of six diffusively coupled Rossler 
systems, (a) The behaviour of xi — X2 when changing the coupling parameter e. This difference is 
chosen as a representative difference of variables associated with the full synchrony subspace A = 
{m = • • • = U6}. The fully synchronous subspace A is stable between e « 0.2 and e « 1.0. (b) The 
behaviour of xi — X3 when changing the coupling parameter e, chosen as a representative difference 
of the variables associated with the synchrony subspace Ax = {ui = 113 = 115,112 = 114 = U6}. 
This partial synchrony subspace is observed above e « 1.0 (when the fully synchrony subspace looses 
stability). The other parameter values are a = 0.2 7 b = 0.2, c = 5.7. For each fixed parameter value 
e (in steps of 0.01 for < e < 1.2), 500 different initial conditions are generated and the difference 
x\ —X2 or x\ — X3 is plotted using the state variables after 40000 iterates with the time step h = 0.01. 
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Table 5.2 

Symbolic adjacency matrix for network in Figure \5.3\ which represents different weights on 
arrows, and the associated coupled cell system. 



Rank 1 




Fig. 5.3. Five globally coupled Lorenz systems and the associated lattice of balanced equivalence 
relations, which has a unique non-trivial balanced equivalence relation (12) (3) (4) (5). 



can induce intermittent bursting for some systems. This bubbling phenomenon is 
observed in synchrony subspaces corresponding to balanced coloring (see the example 
system (14.1) in [18]). 

5.3. Coupled neurons on a random network. The Aldis [4] or Belykh and 
Hasler [9] algorithm finds the minimal balanced coloring, which is the balanced equiv- 
alence relation with the minimal number of colors (i.e. the minimal number of the 
synchronized clusters), and thus the top lattice node. Belykh and Hasler demon- 
strated this using a coupled identical Hindmarsh-Rose model [21 with 30 neurons 
generated by randomly choosing bidirectional identical couplings between any two 
nodes with a small probability. Even though this network has only one type of cell 
and one type of coupling, it is not regular. This network has no apparent symmetry 
using the circular layout (Figure |5.4[ a)) which hides a local reflectional symmetry 
(Figure |5^ b)). The minimal balanced coloring has 23 colors (i.e. 23 synchronized 
clusters). Our algorithm shows the lattice of balanced equivalence relations contains 
only two lattice points, the trivial bottom lattice node (all distinct) and this one non- 
trivial balanced coloring. Modifying this network by adding random edges or rewiring 
existing edges can generate a more complex lattice, for example ten lattice nodes with 



a minimal balanced coloring of 18 clusters (Figure 5.4 'c)). The Python code in the 
supplementary materials finds the lattices for both of these networks. A systematic 
exploration of the lattices of random networks, such as those generated by rewiring 
or the Watts- Strogatz model [40], is an area for future work. 

5.4. Excitatory /Inhibitory coupled neurons. The network of inhibitory 



coupled FitzHugh-Nagumo model neurons shown in Figure 5.5 'a) is discussed in 
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(a) 



(c) 



Fig. 5.4. Thirty neuron networks with bidirectional coupling. Colored nodes represent clusters 
in the minimal balanced coloring, with uncolored (white) nodes for distinct uni-clusters. Here (a) 
and (b) show the same network, from Figure 3(b) in fffi. with 23 clusters. The final figure (c) shows 
a rewired network connecting nodes 13 and 27 in place of 4 and 17, resulting in a richer lattice of 
ten nodes, with 18 clusters in the minimal balanced coloring. Removing the coupling between 4 and 
17 increases the local symmetry, as does coupling nodes 13 and 27. 



[30j [31] . Table 5.3 shows the the 9x9 adjacency matrix C = (cij) with integer 



entries and the associated coupled cell system. A very similar neural network topol- 
ogy (deleting the arrow from neuron 1 to 5) is studied in [L3j using a discrete map 
instead of ODEs. We remark that our algorithm shows both network structures have 
the same lattice of 27 balanced equivalence relations (see lattice generated by the 
Python code in the supplementary materials). The top lattice node is the synchro- 
nized cluster pattern (19) (2378) (46) (5) which was discussed in [13]. 





Fig. 5.5. Nine neurons connected with (a) one coupling type as in [30, 3JS, (b) two coupling 
types (excitatory and inhibitory). Solid arrows represent inhibitory coupling, dashed lines excitatory. 
In (b) the arrows from nodes 3 and 8 are excitatory. 



This network structure was studied as a winnerless competition network [30j [31] , 
in which cluster states (unstable saddle states) are connected along a heteroclinic 
orbit. Such cluster states can correspond to balanced polydiagonals (invariant sub- 
spaces) , which are defined by balanced equivalence relations [8] [7] . Thus the lattice 
of balanced equivalence relations might potentially be used to elucidate the possible 
robust heteroclinic cycles. 
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Table 5.3 

Adjacency matrix which represents the network topology of 9 coupled neurons in Figure \5~5]f a), 
and the associated coupled cell system. Note that this network is not regular since it has three 
input equivalence classes. Note also that the single coupling type constrains the linearized external 
couplings of the three maps f, g and h to be the same. 



To demonstrate multiple arrow types, we modified the previous example to con- 
sider two coupling types, excitatory or inhibitory, as in Figure [K5^ b). By changing the 
outputs of neurons 3 and 8 to be excitatory (i.e. changing four couplings), the number 
of balanced equivalence relations decreased from 27 to 15 (see lattices generated by 
the Python code in the supplementary materials). Table [5^4] shows the corresponding 
symbolic adjacency matrix and the associated coupled cell system. 

symbolic adjacency matrix coupled cell system 
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Table 5.4 

Symbolic adjacency matrix which represents the network topology of 9 coupled neurons with 
excitatory (symbol 'a') and inhibitory (symbol e b') coupling in Figure \5.5}f b), and the associated 
coupled cell system. 



6. Conclusions. Networks in real world applications are inhomogeneous. In- 
dividual systems in a network play different roles and they interact with each other 
in various ways. For example even in a simplified representation of gene regulatory 
networks, genes or proteins can interact either by activation or inhibition. We en- 
coded different types of interaction in such inhomogeneous networks using a symbolic 
adjacency matrix. We considered possible partial synchronies of a given network, 
where the network elements can be grouped into clusters whose dynamics are self- 
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synchronous. We are particularly interested in partial synchronies which are solely 
determined by the structure (topology) of the network, rather than the specific dynam- 
ics (such as function forms and parameter values). Such robust patterns of synchrony 
are associated with balanced equivalence relations, which can be determined by a ma- 
trix computation on the symbolic adjacency matrix. These symbolic adjacency matri- 
ces can alternatively be expressed as a linear combination of integer entry adjacency 
matrices for each coupling type, and the matrix computation applied to each arrow 
type matrix individually (as in the provided Python program). The later is simpler 
to implement as most programming languages or software packages do not support 
symbolic matrices. Using the Symbolic Python library (http://www.sympy.org/) our 
example program can be modified to work on symbolic matrices (not shown as the 
extra dependency complicates installation for no practical benefit). 

The symbolic adjacency matrix therefore specifies the set of balanced equivalence 
relations for a network. From these the refinement relation gives a complete lattice. 
Rather than obtaining the lattice in this way by exhaustive computation, for the spe- 
cial case of regular networks with simple eigenvalues, Kamei [22 , 23 , 24 showed how 
to construct the lattice from building blocks related to the eigenvector /eigenvalues of 
the adjacency matrix, and use it to predict the existence of codimension-one steady- 
state bifurcation branches from the fully synchronous state, and to classify synchrony- 
breaking bifurcation behaviors. Results in [16 , 28, lJ[3j[T7] also relate the eigenvalues 
of the Jacobian of a coupled cell system with the eigenvalues of the adjacency matrix 
of a homogeneous network for synchrony-breaking bifurcation analysis. The connec- 
tions between the algebraic properties of the lattice and the network dynamics are 
potentially of wide interest. 

In general however, such theoretical approaches for the explicit construction of 
the lattice do not yet exist, leaving the "brute force" approach of calculating all 
the possible balanced equivalence relations as the only currently viable route. We 
hope that use of this algorithm will facilitate further theoretical work, and stimulate 
investigation linking lattice properties and synchronous dynamics, and ultimately 
links between network structure and dynamics. 

Acknowledgments. H. Kamei thanks Prof. Ian Stewart for his supervision 
while the foundations of this work was carried out at the University of Warwick. 
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7. Appendix - Algorithm to find the top lattice node. The algorithm we 
use is a generalization of Belykh and Hasler [9] to consider multiple arrow types, or 
equivalently a simplification of the Aldis [4] algorithm where phase one is eliminated 
at the cost of counting absent arrow type/tail- node color combinations. This simpli- 
fication makes it easier to implement than the original Aldis algorithm, yet it is still 
more than fast enough for our needs. The nine-cell network in Figure 5.5 'b) with two 
arrow types is used as an example: 

Step 0. Start by assigning the same color (node class) to each node, here shown 
in red. If multiple node types are considered as in Aldis, then each node type would 
be allocated a unique color. Aldis also classifies the arrows but we skip with that. 

Step 1. In each step compute the "input driven refinement" by tallying the inputs 
to each node according to the color of the node the input is from, and the arrow type. 
After tabulation, unique input combinations give the next node partition. 

See Figure [lT] Here we have one node color (red), and two arrow types (solid 
and dashed), so for each node there are two input counts (solid from red, dashed from 
red). Here Aldis would also compute the same two input counts per node. 




Old partition: (123456789) 



Old color: 










© 


© 


© 


© 


© 


© 


Total: 


(Red)— ► 

\ — y solid 


1 


1 


1 


2 


4 


1 


1 


1 




12 


(Red )---■► 

V_y Dash 








1 




2 








4 


New color: 





© 


© 


© 


© 


© 


© 


© 


© 





New partition: (12378) (4) (5) (6) (9) 

Fig. 7.1. Finding Top Lattice Node (Step 1), input driven refinement of partition (123456789). 



We observe five unique input combinations, and so assign them five colors as the 
"input driven refinement" of the (trivial) input partition. For instance, nodes with one 
solid input from a red node only have been assigned the new partition color orange. 

Step 2. There are now five node colors (shown here as orange, blue, green, yellow 
and cyan). Thus with two arrow types (solid and dashed), we consider ten input types 
(5x2 = 10; solid from orange, . . ., solid from cyan, dashed from orange, . . ., dashed 



from cyan). See Figure 7.2 



We observe six unique input combinations, giving six colors in the new node 
partition. The colors shown are arbitrary, and in the implementation are simply 
integers assigned incrementally. For this example we have reused the colors blue, 
green, yellow and cyan since those node groupings are unchanged. The former orange 
nodes have now been divided into pink and purple nodes. 

Notice that of the ten possible input types tabulated here, the last four are absent. 
The Aldis algorithm avoids counting these. 

Step 3. There are now six node colors (pink, purple, blue, green, yellow and 
cyan), so with two arrow types (solid and dashed) we consider twelve input types 
(6x2 = 12; solid from pink, . . ., solid from cyan, dashed from pink, . . ., dashed from 
cyan). See Figure 7.3 At this iteration the partition of nodes is unchanged, and the 



algorithm halts. This gives the top lattice node. 
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Finding Top Lattice Node (Step 2), refinement of partition (12378) (4) (5) (6) (9) . 
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. Finding Top Lattice Node (Step 3), halts at partition (1)(2378)(4)(5)(6)(9). 
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As in the previous step, some of the possible input types tabulated do not occur 
(five out of twelve), and the Aldis algorithm avoids counting these. 

Remark. When comparing the tables in steps 2 and 3, reusing the same color 
for node groups preserved between iterations highlights that the blue, green, yellow 
and cyan rows are unchanged. In this example at step 3 only the new pink and purple 
rows need be calculated (replacing the orange rows in step 2). This suggests a possible 
speed optimization when finding the top lattice node. 

Differences between Aldis (2008) and our implementation. In the above 
we have noted that the Aldis algorithm avoids computing the zero rows present in 
our tally tables. This is done by an additional phase in each iteration which tracks 
the arrow type and tail node color combinations as arrow equivalence classes, shown 



schematically in Figure 7.4 This figure shows ancestry trees of the node partitions 
(left) and the observed combinations of arrow type (solid or dashed) with tail node 
color (right). By tracking the observed arrow type and tail- node color combinations 
explicitly, absent potential combinations need not be counted (i.e. dashed arrows 
from blue, green, yellow, cyan or pink nodes). 
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Fig. 7.4. Tree of node and arrow partitions as used in the Aldis algorithm. The node partitions 
on the left are (123456789) at step 1 (all red), (12378) (4) (5) (6) (9) at step two (five colors), and finally 
(1) (2378) (4) (5) (6) (9) at step 3 (six colors). The arrow partitions are shown in the disjoint tree on 
the right, using our notation combining the arrow type (solid or dashed) and tail-node partition 
color. The numbers in each box indicate the number of arrows of that type from nodes of that color, 
for example in step 1 there are 12 solid arrows from red nodes. 



The solid arrow tree and the dashed arrow trees in Figure 7.4 (right) are both sub- 
trees of the node partition tree (left). Our approach can be viewed as implicitly using 
the full node partition tree for each arrow type, at the cost of including redundant zero 
branches. This is a tradeoff between algorithmic complexity (Aldis) versus additional 
memory and computational overhead (not noticeable on the graph sizes considered). 
We expect the approach of Aldis to be most beneficial in large networks with many 
arrow types. 



